# This file estimates effects of having ones first two children be of same sex for mothers

# install.packages("dplyr")
# install.packages("date")
# install.packages("AER")

library(dplyr)
library(date)
library(AER)

# Load data
load("mothers_firsttwo.rdata")

# Filter of if age of oldest is unkwown (select parents )
mothers <-
  mothers %>%
  filter(!is.na(age_oldest)) 

# split data file to each election 
# 2009 election
# Subset and create vote variable
# filter of parents with less than two children at the election and a firstborn older than 26

mothers09 <- 
  mothers %>%
  filter(! is.na(stemte_2009) ) %>%
  filter(age_oldest >= 5113 & no_children_09 > 1) %>%
  mutate(vote = stemte_2009 == "Stemte",
         no_children = no_children_09)

#create age dummies for mothers 
#dichotomize twin first
mothers09 <-
  within(mothers09, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })


mothers09 <-
  mothers09 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))


## Repeat for 13
mothers13 <- 
  mothers %>%
  filter(! is.na(stemt) ) %>%
  filter(age_oldest >= 6574 & no_children_13 > 1) %>%
  mutate(vote = stemt == "Stemte",
         no_children = no_children_13)

mothers13 <-
  within(mothers13, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })

mothers13 <-
  mothers13 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))

## Repeat for 14
mothers14 <- 
  mothers %>%
  filter(! is.na(ep_stemt) ) %>%
  filter(age_oldest >= 6574 & no_children_14 > 1) %>%
  mutate(vote = ep_stemt,
         no_children = no_children_14)

mothers14 <-
  within(mothers14, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })

mothers14 <-
  mothers14 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))

# Run models 
# Simple models. 
mod09_r <- lm(vote ~ first_two_same, data = mothers09)
mod13_r <- lm(vote ~ first_two_same, data = mothers13)
mod14_r <- lm(vote ~ first_two_same, data = mothers14)

# Models with fixed effect for mothers and children's age
mod09 <- lm(vote ~ first_two_same + as.factor(age) * as.factor(child_age), data = mothers09)
mod13 <- lm(vote ~ first_two_same + as.factor(age) * as.factor(child_age), data = mothers13)
mod14 <- lm(vote ~ first_two_same + as.factor(age) * as.factor(child_age), data = mothers14)

## Instrument number of children with first child being twin
mod09_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     first_two_same + as.factor(age) * as.factor(child_age), data = mothers09)
mod13_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     first_two_same + as.factor(age) * as.factor(child_age), data = mothers13)
mod14_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     first_two_same + as.factor(age) * as.factor(child_age), data = mothers14)

## Run first-stage models 
mod09_fs <- lm(no_children ~ first_two_same + as.factor(age) * as.factor(child_age), data = mothers09)
mod13_fs <- lm(no_children ~ first_two_same + as.factor(age) * as.factor(child_age), data = mothers13)
mod14_fs <- lm(no_children ~ first_two_same + as.factor(age) * as.factor(child_age), data = mothers14)


# Compile results 

table_mom_09 <-
  rbind(summary(mod09_r)$coefficients,
        summary(mod09)$coefficients[2,],
        summary(mod09_iv)$coefficients[2,])

table_mom_13 <-
  rbind(summary(mod13_r)$coefficients,
        summary(mod13)$coefficients[2,],
        summary(mod13_iv)$coefficients[2,])

table_mom_14 <-
  rbind(summary(mod14_r)$coefficients,
        summary(mod14)$coefficients[2,],
        summary(mod14_iv)$coefficients[2,])

tab_first_stage <- 
  rbind(summary(mod09_fs)$coefficients[1:2,],
        summary(mod13_fs)$coefficients[1:2,],
        summary(mod14_fs)$coefficients[1:2,])

outtab_mom <- 
  rbind(cbind(table_mom_09, with(mothers09, table(first_two_same))),
        cbind(table_mom_13, with(mothers13, table(first_two_same))),
        cbind(table_mom_14, with(mothers14, table(first_two_same))),
        cbind(tab_first_stage, with(mothers14, table(first_two_same))))

write.table(outtab_mom, file = "mother_firsttwo_res.txt")

# Find F diagnostics 
sum_iv09 <- summary(mod09_iv, diagnostics = TRUE)
sum_iv13 <- summary(mod13_iv, diagnostics = TRUE)
sum_iv14 <- summary(mod14_iv, diagnostics = TRUE)

f_diag_mothers <- 
  rbind(sum_iv09$diagnostics[1,],
        sum_iv13$diagnostics[1,],
        sum_iv14$diagnostics[1,])

write.table(f_diag_mothers, file = "f_stat_mothers_firsttwo.txt")

